<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of mvk2</title>
  <meta name="keywords" content="mvk2">
  <meta name="description" content="Two-sided DMRG fast matrix-by-vector product">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../index.html">Home</a> &gt;  <a href="#">tt2</a> &gt; <a href="index.html">@tt_matrix</a> &gt; mvk2.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../index.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for tt2/@tt_matrix&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>mvk2
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>Two-sided DMRG fast matrix-by-vector product</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [y,swp]=mvk2(a,x,eps,nswp,z,rmax,varargin) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">Two-sided DMRG fast matrix-by-vector product
   [Y,SWP]=MVK2(A,X,EPS,[NSWP],[Y],[RMAX],[OPTIONS]) Two-sided DMRG (mvk
   is one-sided). Matrix-by-vector product of a TT-matrix A
   by a TT-tensor X with accuracy EPS. Also, one can specify the number of
   sweeps NSWP, initial approximation Z and the maximal TT-rank RMAX (if
   they become too large)


 TT-Toolbox 2.2, 2009-2012

This is TT Toolbox, written by Ivan Oseledets et al.
Institute of Numerical Mathematics, Moscow, Russia
webpage: http://spring.inm.ras.ru/osel

For all questions, bugs and suggestions please mail
ivan.oseledets@gmail.com
---------------------------</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="../../tt2/@qtt_tucker/diag.html" class="code" title="function [qt]=diag(qt)">diag</a>	Diagonal of a matrix or diagonal matrix from a vector in QTT-Tucker</li><li><a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the QTT-Tucker</li><li><a href="../../tt2/@qtt_tucker/rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>	Computes ranks of the QTT-Tucker tensor</li><li><a href="../../tt2/@qtt_tucker/round.html" class="code" title="function [tt]=round(tt,varargin)">round</a>	Approximate QTT-Tucker with another one with specified accuracy</li><li><a href="conj.html" class="code" title="function [b]=conj(a)">conj</a>	Complex conjugate of a TT-matrix</li><li><a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>	Extract the diagonal of the TT-matrix</li><li><a href="mvk2.html" class="code" title="function [y,swp]=mvk2(a,x,eps,nswp,z,rmax,varargin)">mvk2</a>	Two-sided DMRG fast matrix-by-vector product</li><li><a href="ndims.html" class="code" title="function d = ndims(tt)">ndims</a>	Number of dimensions of a TT-matrix</li><li><a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>	Matrix norm of the TT-matrix</li><li><a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>	All or individual TT-ranks of a TT-matrix</li><li><a href="round.html" class="code" title="function [tt]=round(tt,eps,rmax)">round</a>	Approximate TT-matrix with another one with specified accuracy</li><li><a href="size.html" class="code" title="function [sz] = size(tt)">size</a>	Mode sizes of the TT-matrix</li><li><a href="../../tt2/@tt_tensor/conj.html" class="code" title="function [tt1]=conj(tt)">conj</a>	Compute a complex conjugate of TT-tensor</li><li><a href="../../tt2/@tt_tensor/diag.html" class="code" title="function [tm]=diag(tt)">diag</a>	Constructs diagonal TT-matrix from TT-tensor</li><li><a href="../../tt2/@tt_tensor/ndims.html" class="code" title="function d = ndims(tt)">ndims</a>	Number of dimensions of a tensor</li><li><a href="../../tt2/@tt_tensor/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>	Left and right orthogonalization of the TT-format</li><li><a href="../../tt2/@tt_tensor/rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>	Returns TT-ranks of the given tensor in TT-format</li><li><a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>	Reshape of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/round.html" class="code" title="function [tt]=round(tt,varargin)">round</a>	Approximate TT-tensor with another one with specified accuracy</li><li><a href="../../tt2/@tt_tensor/size.html" class="code" title="function [sz] = size(tt,dim)">size</a>	Mode sizes of the TT-tensor</li><li><a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>	Truncation by absolution precision in Frobenius norm</li><li><a href="../../tt2/core/reort.html" class="code" title="function [y]=reort(u,uadd)">reort</a>	Golub-Kahan reorthogonalization</li><li><a href="../../tt2/core/tt_rand.html" class="code" title="function [tt]=tt_rand(n,d,r)">tt_rand</a>	Generates a random tensor</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="mvk2.html" class="code" title="function [y,swp]=mvk2(a,x,eps,nswp,z,rmax,varargin)">mvk2</a>	Two-sided DMRG fast matrix-by-vector product</li></ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [y,swp]=mvk2(a,x,eps,nswp,z,rmax,varargin)</a>
0002 <span class="comment">%Two-sided DMRG fast matrix-by-vector product</span>
0003 <span class="comment">%   [Y,SWP]=MVK2(A,X,EPS,[NSWP],[Y],[RMAX],[OPTIONS]) Two-sided DMRG (mvk</span>
0004 <span class="comment">%   is one-sided). Matrix-by-vector product of a TT-matrix A</span>
0005 <span class="comment">%   by a TT-tensor X with accuracy EPS. Also, one can specify the number of</span>
0006 <span class="comment">%   sweeps NSWP, initial approximation Z and the maximal TT-rank RMAX (if</span>
0007 <span class="comment">%   they become too large)</span>
0008 <span class="comment">%</span>
0009 <span class="comment">%</span>
0010 <span class="comment">% TT-Toolbox 2.2, 2009-2012</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%This is TT Toolbox, written by Ivan Oseledets et al.</span>
0013 <span class="comment">%Institute of Numerical Mathematics, Moscow, Russia</span>
0014 <span class="comment">%webpage: http://spring.inm.ras.ru/osel</span>
0015 <span class="comment">%</span>
0016 <span class="comment">%For all questions, bugs and suggestions please mail</span>
0017 <span class="comment">%ivan.oseledets@gmail.com</span>
0018 <span class="comment">%---------------------------</span>
0019 n=a.n;
0020 m=a.m;
0021 att=a.tt;
0022 corea=att.core;
0023 psa=att.ps;
0024 ra=att.r;
0025 d=att.d;
0026 <span class="comment">%start_val='fast_mv';</span>
0027 start_val=<span class="string">'rough_mv'</span>;
0028 <span class="comment">%start_val='random';</span>
0029 rmax_loc=8; <span class="comment">%For the &quot;rough&quot; matrix-by-vector product</span>
0030 <span class="keyword">if</span> ( nargin &lt;= 5 || isempty(rmax) )
0031    rmax=1000;
0032 <span class="keyword">end</span>
0033 <span class="keyword">if</span> ( nargin &lt;= 4 || isempty(z) )
0034     
0035     <span class="keyword">if</span> ( strcmp(start_val,<span class="string">'random'</span>) )
0036     rz1=<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(a,1)*<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(x,1);
0037     rzd=<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(a,d+1)*<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(x,d+1);
0038     kf=5;
0039     rz=[rz1;kf*ones(d-1,1);rzd];
0040     z=<a href="../../tt2/core/tt_rand.html" class="code" title="function [tt]=tt_rand(n,d,r)">tt_rand</a>(n,<a href="ndims.html" class="code" title="function d = ndims(tt)">ndims</a>(x),rz);
0041     <span class="keyword">elseif</span> (strcmp(start_val,<span class="string">'rough_mv'</span>))
0042         rmax_loc=8; <span class="comment">%</span>
0043         xloc=<a href="round.html" class="code" title="function [tt]=round(tt,eps,rmax)">round</a>(x,0,rmax_loc);
0044         aloc=<a href="round.html" class="code" title="function [tt]=round(tt,eps,rmax)">round</a>(a,0,rmax_loc);
0045         z=<a href="round.html" class="code" title="function [tt]=round(tt,eps,rmax)">round</a>(aloc*xloc,eps); 
0046         <span class="comment">%z=aloc*xloc;</span>
0047     <span class="keyword">elseif</span> (strcmp(start_val,<span class="string">'fast_mv'</span>) )
0048          rz1=<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(a,1)*<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(x,1);
0049       rzd=<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(a,d+1)*<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(x,d+1);
0050       kf=5;
0051       rz=[rz1;kf*ones(d-1,1);rzd];
0052       z=<a href="../../tt2/core/tt_rand.html" class="code" title="function [tt]=tt_rand(n,d,r)">tt_rand</a>(n,<a href="ndims.html" class="code" title="function d = ndims(tt)">ndims</a>(x),rz);
0053        z=<a href="mvk2.html" class="code" title="function [y,swp]=mvk2(a,x,eps,nswp,z,rmax,varargin)">mvk2</a>(a,x,max(eps,1e-2),10,z,rmax); <span class="comment">%First, do it with bad accuracy</span>
0054     <span class="keyword">end</span>
0055 <span class="keyword">end</span>
0056 <span class="keyword">if</span> ( nargin &lt;= 3 || isempty(nswp) )
0057   nswp = 40;
0058 <span class="keyword">end</span>
0059 y=z;
0060 
0061 <span class="comment">%Parameters section</span>
0062 kick_rank=6;
0063 verb=false;
0064 
0065 <span class="comment">%Warmup is to orthogonalize Y from right-to-left and compute psi-matrices</span>
0066 <span class="comment">%for Ax</span>
0067 psi=cell(d+1,1); <span class="comment">%Psi-matrices</span>
0068 <span class="comment">%psi{d+1}=1; psi{1}=1;</span>
0069 psi{d+1}=eye(<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(y,d+1)); 
0070 psi{1}=eye(<a href="rank.html" class="code" title="function [r]=rank(a,varargin)">rank</a>(y,1));
0071 <span class="comment">%Here we will add convergence test</span>
0072 <span class="comment">%Warmup: right-to-left QR + computation of psi matrices</span>
0073 
0074 corex=x.core;
0075 psx=x.ps;
0076 rx=x.r;
0077 
0078 swp=1;
0079 converged=false;
0080 <span class="keyword">while</span> (swp &lt;= nswp &amp;&amp; ~converged)  
0081 psy=y.ps;
0082   ry=y.r;
0083   corey=y.core;
0084      pos1=psy(d+1);
0085 cr1=corey(psy(d):psy(d+1)-1);
0086 
0087 <span class="keyword">for</span> i=d:-1:2  
0088    cr2=corey(psy(i-1):psy(i)-1);
0089    cr1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[ry(i),n(i)*ry(i+1)]);
0090    cr2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[ry(i-1)*n(i-1),ry(i)]);
0091    cr1=cr1.';
0092    [q,rm]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(cr1,0); rn=<a href="size.html" class="code" title="function [sz] = size(tt)">size</a>(q,2); rm=rm.';
0093    q=q.'; 
0094    ry(i)=rn;
0095    corey(pos1-ry(i+1)*n(i)*ry(i):pos1-1)=q(:);
0096    <span class="comment">%Convolution is now performed for psi(i) using psi(i+1) and corea, and</span>
0097    <span class="comment">%(new) core q</span>
0098    cra=corea(psa(i):psa(i+1)-1); cra=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra,[ra(i),n(i),m(i),ra(i+1)]);
0099    cry=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(<a href="conj.html" class="code" title="function [b]=conj(a)">conj</a>(q),[ry(i),n(i),ry(i+1)]);
0100    crx=corex(psx(i):psx(i+1)-1); crx=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx,[rx(i),m(i),rx(i+1)]);
0101    pscur=psi{i+1}; pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[ra(i+1),rx(i+1),ry(i+1)]); <span class="comment">%ra,rx,ry</span>
0102    <span class="comment">%First, convolve over rx(i+1)</span>
0103    crx=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx,[rx(i)*m(i),rx(i+1)]);
0104    pscur=permute(pscur,[2,1,3]); pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i+1),ra(i+1)*ry(i+1)]);
0105    pscur=crx*pscur; <span class="comment">%pscur is now rx(i)*m(i)*ra(i+1)*ry(i+1)</span>
0106    <span class="comment">%Convolve over m(i),ra(i+1),n(i),ry(i+1)</span>
0107     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i),m(i)*ra(i+1),ry(i+1)]);
0108     pscur=permute(pscur,[1,3,2]); 
0109     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i)*ry(i+1),m(i)*ra(i+1)]);
0110     cra=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra,[ra(i)*n(i),m(i)*ra(i+1)]); cra=cra.';  
0111     pscur=pscur*cra; 
0112     <span class="comment">%pscur is now rx(i)*ry(i+1)*ra(i)*n(i), it is left to convolve over</span>
0113     <span class="comment">%n(i)*ry(i+1)</span>
0114     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i),ry(i+1),ra(i),n(i)]);
0115     pscur=permute(pscur,[3,1,4,2]);
0116     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i)*ra(i),n(i)*ry(i+1)]);
0117     cry=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry,[ry(i),n(i)*ry(i+1)]); cry=cry.';
0118     pscur=pscur*cry;
0119     psi{i}=pscur;
0120    <span class="comment">%End of psi-block</span>
0121    pos1=pos1-ry(i+1)*n(i)*ry(i);
0122    cr1=cr2*rm;
0123    
0124 <span class="keyword">end</span>
0125 corey(pos1-ry(2)*n(1)*ry(1):pos1-1)=cr1(:);
0126 pos1=pos1-ry(2)*n(1)*ry(1);
0127 corey=corey(pos1:numel(corey)); <span class="comment">%Truncate unused elements</span>
0128 
0129   
0130 
0131   <span class="comment">%left-to-right dmrg sweep</span>
0132   pos1=1;
0133   cry_old=corey;
0134   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0135   converged=true;
0136   ermax=0;
0137   <span class="keyword">for</span> i=1:d-1
0138      <span class="comment">%We care for two cores, with number i &amp; number i+1, and use</span>
0139      <span class="comment">%psi(i) and psi(i+2) as a basis; also we will need to recompute</span>
0140      <span class="comment">%psi(i+1)</span>
0141      ps1=psi{i}; ps2=psi{i+2};
0142      cra1=corea(psa(i):psa(i+1)-1); cra2=corea(psa(i+1):psa(i+2)-1);
0143      crx1=corex(psx(i):psx(i+1)-1); crx2=corex(psx(i+1):psx(i+2)-1);
0144      <span class="comment">%our convolution is</span>
0145      <span class="comment">%ps1(ra(i),rx(i),ry(i))*cra1(ra(i),n(i),m(i),ra(i+1))*</span>
0146      <span class="comment">%cra2(ra(i+1),n(i+1),m(i+1),ra(i+2))*</span>
0147      <span class="comment">%*ps2(ra(i+2),rx(i+2),ry(i+2))</span>
0148      <span class="comment">%*crx1(rx(i),m(i),rx(i+1))*cr2x(rx(i+1)*m(i+1)*rx(i+2))</span>
0149      <span class="comment">%Scheme ps1*crx1 over rx(i)</span>
0150      
0151      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i),rx(i),ry(i)]); 
0152      ps1=permute(ps1,[1,3,2]); ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i)*ry(i),rx(i)]);
0153      crx1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx1,[rx(i),m(i)*rx(i+1)]);
0154      ps1=ps1*crx1; <span class="comment">%ps1 is now ra(i)*ry(i)*m(i)*rx(i+1)</span>
0155      <span class="comment">%Now convolve with matrix A over ra(i)*m(i)</span>
0156      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i),ry(i),m(i),rx(i+1)]);
0157      ps1=permute(ps1,[2,4,1,3]);
0158      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i)*rx(i+1),ra(i)*m(i)]);
0159      cra1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra1,[ra(i),n(i),m(i),ra(i+1)]);
0160      cra1=permute(cra1,[1,3,2,4]); 
0161      cra1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra1,[ra(i)*m(i),n(i)*ra(i+1)]);
0162      ps1=ps1*cra1; <span class="comment">%ps1 is now ry(i)*rx(i+1)*n(i)*ra(i+1)</span>
0163      <span class="comment">%Then the ``same'' convolution is carried over for second pair of</span>
0164      <span class="comment">%cores</span>
0165      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[ra(i+2),rx(i+2),ry(i+2)]);
0166      ps2=permute(ps2,[2,1,3]);
0167      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+2),ra(i+2)*ry(i+2)]);
0168      crx2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx2,[rx(i+1)*m(i+1),rx(i+2)]);
0169      ps2=crx2*ps2; <span class="comment">%ps2 is now rx*(i+1)*m(i+1)*ra(i+2)*ry(i+2)</span>
0170      <span class="comment">%Convolve over m(i+1)*ra(i+2)</span>
0171      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1),m(i+1)*ra(i+2),ry(i+2)]);
0172      ps2=permute(ps2,[2,1,3]);
0173      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[m(i+1)*ra(i+2),rx(i+1)*ry(i+2)]);
0174      cra2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra2,[ra(i+1)*n(i+1),m(i+1)*ra(i+2)]);
0175      ps2=cra2*ps2; 
0176      <span class="comment">%ps2 is now ra(i+1)*n(i+1)*rx(i+1)*ry(i+2)</span>
0177      <span class="comment">%Now form superblock by contraction ps1 &amp; ps2</span>
0178      <span class="comment">%over ra(i+1)*rx(i+1)</span>
0179      ps0=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
0180      ps0=permute(ps0,[1,3,2,4]);
0181      ps0=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps0,[ry(i)*n(i),rx(i+1)*ra(i+1)]);
0182      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[ra(i+1),n(i+1),rx(i+1),ry(i+2)]);
0183      ps2=permute(ps2,[3,1,2,4]); 
0184      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1)*ra(i+1),n(i+1)*ry(i+2)]);
0185      super_core=ps0*ps2; <span class="comment">%super_core is ry(i)*n(i)*n(i+1)*ry(i+2)</span>
0186     <span class="keyword">if</span> ( i &gt; 1 )
0187      <span class="comment">%Compute previous supercore</span>
0188         cr1=corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1);
0189         cr2=cry_old(psy(i+1):psy(i+2)-1); 
0190         cr1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[ry(i)*n(i),ry(i+1)]);
0191         cr2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[ry(i+1),n(i+1)*ry(i+2)]);
0192          super_core_old=cr1*cr2; 
0193          er=<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(super_core_old(:)-super_core(:))/<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(super_core(:));
0194      
0195         <span class="keyword">if</span> ( er &gt; eps ) 
0196             converged=false;
0197         <span class="keyword">end</span>
0198         ermax=max(er,ermax);
0199      <span class="comment">%if ( verb )</span>
0200      <span class="comment">%  fprintf('i=%d er=%3.2e \n',i,er);</span>
0201      <span class="comment">%end</span>
0202      <span class="keyword">end</span>
0203      [u,s,v]=svd(super_core,<span class="string">'econ'</span>);
0204      s=<a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>(s); 
0205      r=<a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>(s,eps/sqrt(d-1)*<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(s)); r=min(r,rmax);
0206      u=u(:,1:r); s=s(1:r); v=v(:,1:r); v=v*<a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>(s);
0207      
0208      <span class="comment">%Kick rank</span>
0209      
0210      ur=randn(<a href="size.html" class="code" title="function [sz] = size(tt)">size</a>(u,1),kick_rank);
0211      <span class="comment">%Orthogonalize ur to u by Golub-Kahan reorth</span>
0212      u=<a href="../../tt2/core/reort.html" class="code" title="function [y]=reort(u,uadd)">reort</a>(u,ur);
0213      radd=<a href="size.html" class="code" title="function [sz] = size(tt)">size</a>(u,2)-r; 
0214      <span class="keyword">if</span> ( radd &gt; 0 )
0215         vr=zeros(<a href="size.html" class="code" title="function [sz] = size(tt)">size</a>(v,1),radd);
0216         v=[v,vr];
0217      <span class="keyword">end</span>
0218      r=<a href="size.html" class="code" title="function [sz] = size(tt)">size</a>(u,2);
0219      
0220      
0221      ry(i+1)=r;
0222      <span class="comment">%u is ry(i)*n(i)*ry(i+1)</span>
0223      <span class="comment">%core_new(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);</span>
0224      corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:); 
0225      
0226      u=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(u,[ry(i),n(i),ry(i+1)]);
0227 
0228      <span class="comment">%Compute new psi</span>
0229      <span class="comment">%ps1 is ry(i)*rx(i+1)*n(i)*ra(i+1) with u over ry(i)*n(i)</span>
0230      u=<a href="conj.html" class="code" title="function [b]=conj(a)">conj</a>(u);
0231      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
0232      ps1=permute(ps1,[4,2,1,3]);
0233      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i+1)*rx(i+1),ry(i)*n(i)]);
0234      u=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(u,[ry(i)*n(i),ry(i+1)]);
0235      ps1=ps1*u;
0236      psi{i+1}=ps1;
0237      <span class="comment">%Compute (?) new v</span>
0238      pos1=pos1+ry(i)*n(i)*ry(i+1);
0239      v=v';
0240      <span class="comment">%v=reshape(v,[ry(i+1),n(i+1),ry(i+2)]);</span>
0241      corey(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:);
0242   <span class="keyword">end</span>
0243   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0244 y.core=corey;
0245 y.r=ry;
0246 y.ps=psy;
0247 swp=swp+1;
0248 <span class="keyword">if</span> ( verb )
0249 fprintf(<span class="string">'swp=%d er=%3.2e trunk=%3.2e \n'</span>,swp,ermax,eps/sqrt(d-1));
0250 <span class="keyword">end</span>
0251 <span class="keyword">end</span>
0252 <span class="keyword">if</span> ( swp == nswp ) 
0253   fprintf(<span class="string">'mvk2 warning: error is not fixed for maximal number of sweeps %d\n'</span>, swp); 
0254 <span class="keyword">end</span><span class="comment">%end</span>
0255 <span class="comment">%p1=tt_mvdot(core(a),core(x),y0);</span>
0256 <span class="comment">%p2=dot(y,tt_tensor(y0));</span>
0257 <span class="comment">%abs(p1-p2)/abs(p1)</span>
0258 <span class="comment">%keyboard;</span>
0259 
0260 
0261</pre></div>
<hr><address>Generated on Wed 08-Feb-2012 18:20:24 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>